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Abstract 

There is great potential to be explored regarding the use of agent-based modelling and simulation as an alternative 
paradigm to investigate early-stage cancer interactions with the immune system. It does not suffer from some limitations of 
ordinary differential equation models, such as the lack of stochasticity, representation of individual behaviours rather than 
aggregates and individual memory. In this paper we investigate the potential contribution of agent-based modelling and 
simulation when contrasted with stochastic versions of ODE models using early-stage cancer examples. We seek answers to 
the following questions: (1) Does this new stochastic formulation produce similar results to the agent-based version? (2) Can 
these methods be used interchangeably? (3) Do agent-based models outcomes reveal any benefit when compared to the 
Gillespie results? To answer these research questions we investigate three well-established mathematical models describing 
interactions between tumour cells and immune elements. These case studies were re-conceptualised under an agent-based 
perspective and also converted to the Gillespie algorithm formulation. Our interest in this work, therefore, is to establish a 
methodological discussion regarding the usability of different simulation approaches, rather than provide further biological 
insights into the investigated case studies. Our results show that it is possible to obtain equivalent models that implement 
the same mechanisms; however, the incapacity of the Gillespie algorithm to retain individual memory of past events affects 
the similarity of some results. Furthermore, the emergent behaviour of ABIVIS produces extra patters of behaviour in the 
system, which was not obtained by the Gillespie algorithm. 
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Introduction 

In previous work, three case studies using established mathe- 
matical models of immune interactions with early-stage cancer 
were considered in order to investigate the additional contribution 
of ABMS to ODE models simulation [1]. These case studies were 
re-conceptuahsed under an agent-based perspective and the 
simulation results were compared with those from the ODE 
models. Our results showed that, apart from the well known 
differences between these approaches (as those outlined for 
example, in Schieritz and Milling [2]), further insight from using 
ABMS was obtained, such as extra population patterns of 
behaviotu'. 

In this work we apply the Gillespie algorithm [3,4], which is a 
variation of the Monte Carlo method, to create stochastic versions 
of the original ODE models investigated in [1]. We aim to 
reproduce the variability embedded in the ABMS systems to the 
mathematical formulation and verify whether results resemble 
each other. In addition, due to the fact that the Gillespie algorithm 
also regards integer quantities for their elements, we hope that this 
method overcomes some differences observed when comparing 
atomic agents represented in the ABMS with possible fractions of 
elements observable in the ODE results. 

To the best of our knowledge, current literature regarding the 
direct comparison of the Gillespie algorithm and ABMS is scarce. 



We want therefore to answer research questions such as: (1) Does 
this new stochastic formulation produce similar results to ABMS? 
(2) Can these methods be used interchangeably for our case 
studies? (3) Does the stochastic model implemented using the 
Gillespie algorithm also find the extra patterns revealed by the 
ABMS? We aim to establish a methodological discussion regarding 
the benefits of each approach for biological simulation, rather than 
provide further insights into the biological aspects of the problems 
studied. We therefore intend to compare the dynamics of each 
approach and observe the outcomes produced over time. We hope 
that this study provides further insights into the potential usability 
and contribution of ABMS to systems simulation. 

Case Studies 

The case studies used for our comparison regard models with 
different population sizes, varying modelling effort and model 
complexity. We hope that by tackling problems with different 
characteristics, a more robust analysis of our experiments is 
performed. The features of each case study are shown in Table 1. 

The first case study considered is based on an ODE model 
involving interactions between generic tumour and effector cells. 
The second case study adds to the previous model the influence of 
the IL-2 cytokine molecules in the immune responses. The third 
case study comprises a model of interactions between effector cells. 
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Table 1. Case studies considered. 





Case Study 


Number of populations 


Population size 


Complexity 


1) Tumour/Effector 


2 


5 to 600 


Low 


2) Tumour/Efector/IL-2 


3 


10* 


Medium 


3) Tumour/Effector/IL-2/TGF-i? 


4 


10" 


Higii 



doi:10.1371/journal.pone.0095150.t001 



tumour cells, and IL-2 and TGF-/? molecules. For all case studies, 
three approaches are presented: the original mathematical model, 
its conversion into the Gillespie algorithm model and the ABMS 
model. To answer our research questions, the Gillespie and ABMS 
approaches outcomes are compared. Our results show that for 
most cases the Gillespie algorithm does not produce outcomes 
statistically similar to the ABMS. In addition, Gillespie is incapable 
to reproduce the extreme patterns observed in the ABMS 
outcomes for the last case study. 

The remainder of this paper is organized as follows. The next 
section introduces the literature review comparing stochastic ODE 
models and ABMS for different simulation domains. First, we 
show general work that has been carried out in areas such as 
economics and operations research, and then we focus on research 
concerned with the comparison for immunological problems. 
Finally, we discuss gaps in the literature regarding cancer research. 
In the following section we introduce our Gillespie and agent- 
based modelling development processes and the methods used for 
conducting the experimentation. Subsequently we present our case 
studies, comparison results and discussions. In the fmal section we 
draw our overall conclusions and outline future research 
opportunities. 

Related Work 

Current in-silico approaches used in early-stage cancer research 
include computational simulation of compartmental models, 
individual-based models and rule-based models. Compartmental 
models adopt an aggregate representation of the elements in the 
system. They include deterministic methods such as ordinary 
differential equation (ODE) models, system dynamics (SD) models 
and partial differential equation (PDF) models. These models have 
been largely employed in the study of dynamics between cancer 
cells and tumour cells [5,6], therapies for cancer [7], tumour 
responses to low levels of nutrients [8-1 1] and tumour vascular- 
ization[12,13]. Although these models have been very useful to 
understand and uncover various phenomena, they present several 



limitations. For instance, they do not encompass emergent 
behaviour and stochasticity. In addition, it is difficult to keep a 
record of individual behaviour and memory over the simulation 
course [14,15]. Stochastic compartmental models include Monte 
Carlo simulation models, which are computational algorithms that 
perform random sampling to obtain numerical results [16]. 
Amongst others, they are useful for simulating biological systems, 
such as cellular interactions and the dynamics of infectious diseases 
[17]. As these methods rely on stochastic process to produce their 
outputs, they overcome some of the limitations of the deterministic 
compartmental models, as they allow for variability of outcomes. 
The individuals in these models, however, do not have any sort of 
memory of past events. Rule-based models are a relatively new 
research area mostly focused on modelling and simulating 
biochemical reactions, molecular interactions and cellular signal- 
ling. The literature regarding the application of rule-based models 
to interactions between the immune system and cancer cells, 
however, is scarce. Individual-based models, or agent-based 
modelling and simulation (ABMS), relax the aggregation assump- 
tions present in compartmental approaches and allow for the 
observation of the behavior of the single cells or molecules 
involved in the system. This approach has also been applied to 
early-stage cancer research [1,18]. 

The differences between deterministic compartmental models 
and individual-based models are well known in operations 
research [2,19-21] and have also been studied in epidemiology 
[22,23] and system's biology [24-26]. Deterministic compartmen- 
tal models assume continuous values for the individuals in the 
system, whereas in ABMS individual agents are represented. This 
peculiarity of each approach highly impacts the simulation results 
similarity depending on the size of the populations [1]. There is 
still however the need for further investigations between the 
interchangeable use of some Monte Carlo methods and ABMS. 

Approaches Comparison 

As mentioned previously, there are few studies that compare the 
Gillespie algorithm with ABMS. Most of these studies regard 



Table 2. Reactions for case 1. 





Phenomenon 


Reaction equation 


Rate law (per cell) 


Tumour cell birth 


T^lxT 


aT 


Tumour cell death 


2*T^T 


abT^ 


Tumour cell death by effector cells 


T + E~*E 


nTE 


Effector proliferation 


E^2xE;T 


pTE 
i+T 


Effector death by fighting tumour 


T + E~*T 


mTE 


Effector death 


E^ 


dE 


Effector supply 


-^E 


s 



doi:10.1371/journa!.pone.0095150.t002 
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Table 3. Agents' parameters and behaviours for case 1. 





Agent 


Parameters 


Reactive behaviour 


Proactive behaviour 


Tumour Cell 


a and h 


Dies {with age) 






a and b 




Proliferates 




m 




Damages effector cells 




n 


Dies killed by effector cells 




Effector Cell 


m 


Dies (with age) 






d 


Dies per apoptosis 






p and g 




Proliferates 




s 


Is injected as treatment 





doi:10.1371/journal.pone.0095150.t003 



research in economic models and immunology. To the best of our 
knowledge there is no literature regarding the direct comparison of 
these methods to early-stage interactions between the immune 
system and tumour cells. This section describes relevant researches 
in several areas, which provided further insights into the gaps in 
the current literature and the research questions addressed in this 
paper. 

There are a few attempts of re-conceptualizing agent-based 
models into simpler stochastic models of complex systems in 
economics. For instance, Daniunas et al. [27] start from simple 
models with estabhshed agent-based versions (which they named 
"the model's microscopic version") and try to obtain an equivalent 
macroscopic behavior. They consider microscopic and macro- 
scopic versions of the herding model proposed by Kirman [28] 
and the diffusion of new products, proposed by Bass in [29]. They 
conclude that such simple models are easily replicated in a 
stochastic environment. In addition, the authors state that for the 
economics field, only very general models, such as those studied in 
their article, have well established agent-based versions and can be 
described by stochastic or ordinary differential equations. How- 
ever, as the complexity of the microscopic environment increases, 
it becomes challenging to obtain resembling results with stochastic 
simulations and further developments need to be pursuit. 
Furthermore, the authors debate that the ambiguity present in 
the microscopic description in complex systems is an objective obstacle 
for quantitative modeling and needs further studying. 

Stracquadanio et al. [30] investigate the contributions of ABMS 
and the Gillespie method for immune modelling. The authors, 
however, do not apply both methods to the same problem. 
Instead, for the first approach, they chose to investigate a large- 
scale model involving interactions of immune cells and molecules. 



This model's objective was to simulate the immune elements 
interplay over time. For the Gillespie approach, the authors 
investigate a stochastic viral infection model. The authors point 
out three factors that play a major role in the modeling outcome 
when comparing ABMS and Gillespie; simulation time, model 
precision and accuracy, and model applicability. Regarding time, 
the authors state that stochastic models implemented with the 
Gillespie algorithm are preferred. On the other hand, ABMS 
permits more control over simulation runtime as it keeps record of 
the behavior of each single entity involved in the system. 
Regarding applicability, the authors argue that traditional 
Gillespie methods do not account for spatial information, which 
can be detrimental to the model accuracy given the fact that many 
immune interactions occur within specific spatial regions of the 
simulation environment. 

Karkutla [31] compares two biological simulators: GridCell, 
which is a stochastic tool based on Gillespie's, and his new 
developed ABMSim, which is a simulation tool based on ABMS. 
GridCeU was developed to overcome the issues in traditional 
Gillespie's, as pointed out by Stracquadanio et al. [30]. It is a 
stochastic tool able to tackle non-homogeneity effectively by 
addressing issues of crowding and localization. In GridCell, 
however, the problem of tracking individual behaviour and 
determining particular characteristics to each element stiU exists. 
ABMSim was therefore developed to overcome these issues. 
GridCeU has been compared with ABMSim qualitatively and 
quantitatively and the two tools have produced similar results in 
the author experiments. 

In our work we further study the differences between the 
approaches outcomes and investigate whether under different 
problem characteristics for early-stage cancer we still obtain 



Table 4. Transition rates calculations from the mathematical equations for case 1. 



Agent Transition MatKiematlcal equation Transition rate 

Tumour Cell proliferation aTi\ — Tb) a — {TotalTumour.h) 

death aT{\ - Tb) a -(TottilTumour.h) 

dieKilledByEffector nTE n.TotalEffectorCells 

causeEffectorDamage mTE m 

Effector Cell Reproduce p.TouiiTumourCfih 

g+T g+ ToUilTninaurCc'lls 

DiePerAge dE d 

DiePerApoptosis mTE message from tumour 

doi:10.1371/journal.pone.0095150.t004 
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Tumour Cell 



alive 



* ^ causesEffectorPamage 



diekilledByEffector 



death 



proliFeratlon 



^^^dead 

(a) Tumour cell agent 

Figure 1. ABMS state charts for case 1. 

doi:10.1371/journal.pone.0095150.g001 



Effector Cell 



diePerApoptosis 




killTumour 



diePei Age 



(b) Effector cell agent 



similar results. In the next section we introduce the methodology 
used to conduct our investigations. 

Methods 

This section introduces the research methodology used for the 
development of our simulation models and for the experimenta- 
tion performed in the following sections. As mentioned previously, 
our investigations concern the use of three cases to answer 
research questions regarding the application of the Gillespie 
algorithm and ABMS interchangeably for early-stage cancer 
models. For each case study, there is a well established ODE 
model from the literature and its correspondent agent-based 
model, that we previously developed in [1]. The ODE model 
simulations are implemented using the ODE solver module from 
MATLAB (2011) 

The Gillespie algorithm is implemented by the direct conversion 
of the original mathematical equations into reactions and 
simulating them under the COPASI 4.8 {Build'iS) simulator 
environment. The method used for the stochastic simulations is the 
Gillespie algorithm adapted using the next reaction method [32], 
with interval sizes of 0.1, integration interval between 0 and 1 and 
maximum internal steps of 10^. 

ABMS is a modelling and simulation technique that employs 
autonomous agents that interact with each other. The agents' 
behaviour is described by rules that determines how they learn, 
interact and adapt. The overall system behaviour is given by the 
agents individual dynamics as well as their interactions. Our agent- 
based models are implemented using the AnyLogic™ 6.5 
educational version (XJ Technologies 2010) [33]. This approach 
is developed by using state charts and tables containing each agent 
description. The state charts show the different possible states of an 
entity and define the events that cause a transition from one state 



to another. In order to facilitate the understanding of the agent- 
based model, we reproduce here the models developments, which 
were based on [1] (The ABMS and GiUespie models are available 
for download in http://anytips.cs.nott.ac.uk/wiki/index.php/ 
Resources). 

Methodology for Results Comparison 

As the GiUespie and ABMS are both stochastic simulation 
methods, we ran five hundred replications for each case study and 
calculated the mean values for the outputs. For all approaches, the 
rates (for cellular death, birth, etc.) employed were the same as 
those established by the mathematical model. 

In addition, in order to investigate any statically significant 
differences between the ABMS and Gillespie techniques for the 
case studies, we implement a mixed effect model. This is a type of 
regression that considers both fixed and random effects. This 
method accounts for correlation caused by repeating the measure 
over time (i.e., the tumour cell count is correlated over time for 
each simulation run). The mixed effect analysis was implemented 
in the programming langxiage R using the package NLME [34]. 

As a mixed effect model requires finding parameters for a 
regression model, it is not suitable when considering the whole 
time period. This is because in cases 2 and 3 the tumour dynamics 
has a damping oscillation and the function describing this 
dynamics is unknown (see pages 1 1 and 14). Instead, the sequence 
of local maxima and minima are used. It can been seen that these 
are converging and any statistical deviation between these 
sequences for the different simulation techniques indicate differ- 
ences between the output of the techniques. If the simulations from 
the ABMS and the Gillespie technique come from the same 
distribution, then there would be no statistical difference between 
the maxima and minima over time. Therefore, we investigate two 



Table 5. Simulation parameters for the four different scenarios under investigation. 





Scenario 


b 


d 


s 


1 


0.002 


0.1908 


0.318 


2 


0.004 


2 


0.318 


3 


0.002 


0.3743 


0.1181 


4 


0.002 


0.3743 


0 



For the other parameters, the values are the same in all experiments, i.e. ^=1.636, ^ = 20.19, h; = 0. 00311, = ! and p = lA3i. 
doi:l 0.1 371/journal.pone.C0951 50.t005 
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Scenario 1 using the Gillespie Method 
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Figure 2. Results for the first case study. 

doi:1 0.1 371 /journal.pone.00951 50.g002 



null hypotheses. The first is that the function of local maxima is the 
same for the ABMS and the Gillespie algorithm simulations. And 
the second is that the function of local minima is the same for the 
ABMS and the Gillespie algorithm simulations. We use a 1 % 
significance level. 

There is no standard technique for estimating the required 
sample size for non-linear mixed effect models for a defined power 
when the measure of effect is known [35]. Therefore, the 
simulations are run 500 times as this will increase the statistical 
power and increase the probability of a true positive in the 
statistical analysis. A false negative is still possible if there is only a 
small effect size, but if the effect is small, it is of less interest. 

Case 1: Interactions between Tumour Cells and Generic 
Effector Cells 

The first case considers tumour cells growth and their 
interactions with general immune effector cells, as defined in [8]. 
According to the model, effector cells search and kill the tumour 
cells inside the organism. They proliferate proportionally to the 
number of existing tumour cells. As the quantities of effector cells 
increase, their capacity of eliminating tumour cells is augmented. 
Immune cells proliferate and die per apoptosis, which is a 
programmed cellular death. In the model, cancer treatment is also 



considered and it consists of injections of new effector cells in the 
organism. 

Mathematically, the interactions between tumour cells and 
immune effector cells are defmed as follows [8]: 



dT 
dt 

dE 
~di 



■-Tf{T)-dT(T,E) 

--pE(T,E)-dE{T,E) 
-aE(E) + (^{T) 



(1) 



(2) 



where 

• T is the number of tumour cells, 

• £ is the number of effector cells, 

• f (T) is the growth of tumour cells, 

• dT{T,E) is the number of tumour cells killed by effector cells, 

• Pe(T,E) is the proliferation of effector cells, 

• dE(T,E) is the death of effector cells when fighting tumour 



cells. 



Table 6. The fixed parameter values returned by the nnixed-effect model and their significance. 




Method 


Parameter 


Value 


std error 


p-value 


ABS 


a 


0.03393 


0.0004397 


0 


SODE 


a 


0.02925 


0.0006126 


0 


ABS 


b 


41.15847 


0.6149638 


0 


SODE 


h 


52.07957 


1 .0328222 


0 



doi:1 0.1 371 /journal.pone.00951 50.t006 
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Table 7. Reactions for case 2. 





Name 


Reaction equation 


Rate law (per cell) 


Effector cell recruitment 




cT 


Effector cell proliferation 




Si +4) 


Effector cell death 








—*E 


s\ 


Tumour cell birth 


T^lxT 


aT 


Tumour cell death 


2 * T->T 


abT^ 


Tumour cell death by effector cells 


T + E^E-T 


a,TE 
Sl+T 


IL-2 production 


-^I:ET 


es + T 


IL-2 decay 


/-» 


ft 


IL-2 supply 




.s-2 



doi:10.1371/journa!.pone.0095150.t007 



• ueIE) is the death (apoptosis) of efTector cells, 

• <I)(r) is the treatment or influx of cells. aE(E) = dE (7) 

m=s (8) 

Table 2 shows the mathematical equations converted into 
reactions and their respective rate laws per cell. 

In the agent-based model there are two classes of agents, the 
tumor cells and the effector cells, as described in [1]. Table 3 
shows the parameters and behaviours corresponding to each agent 
state. For our agents, state charts are used to represent the 
difiFerent states each entity is in. In addition, transitions are used to 
indicate how the agents move from one state to another. Events 
are also employed and they indicate that certain actions are 
scheduled to occur in the course of the simulation, such as 
injection of treatment. The state chart representing the tumour 
cells is shown in Figure 1(a), in which an agent proliferates, dies 
with age or is killed by effector cells. In addition, tumour cells 



The Kuznetsov model [8] defines the functions f(T), djiT^E), 
Pe(E,T), dE{E,T), ttEiE) and as shown below: 

f(,T) = a(\-bT) (3) 
dT(T,E) = nTE (4) 

PEiE,T)='^ (5) 
dE(E,T) = mTE (6) 



Table 8. Agents' parameters and behaviours for case 2. 





Agent 


Parameters 


Reactive behaviour 


Proactive behaviour 


Effector Cell 


nm2 


Dies 






pi and gl 




Reproduces 




c 


Is recruited 




■vl Is injected as treatment 




pi and g3 




Produces IL-2 




aa and gl 




Kills tumour cells 


Tumour Cell 


a and h 


Dies 






a and h 




Proliferates 




aa and gl 


Dies killed by effector cells 






c 




Induces effector recruitment 


IL-2 


pi and ^3 


Is produced 






muZ 


Is lost 






si 


Is injected 





doi:10.1371/journa!.pone.0095150.t008 
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Table 9. Transition rates calculations from the mathematical equations for case 2. 



Agent 


Transition 


Mathematical equation 


Transition rate 


Effector Cell 


Reproduce 


Pi-IlE 


pi .TotallLJl.TotalEjfector 


gl+ILJl 


g\ + TolallLJl 




Die 








killTumour 


a^ET 
gl+T 


TolalTiimoiir 
g2 + TolaiTumour 




ProducelL2 


plET 


p2. ToliilTiinwiir 
g3 + ToUdTumour 


Tumour Cell 


Reproduce 


aT{\~hT) 


a — iTotalTimiour.h) 




Die 


aT{\-hT) 


a — {TotalTumour .h) 




DieKilledByEffector 


a^TE 
g2+T 


message from effector 


IL-2 


Loss 







doi:10.1371/journal.pone.0095150.t009 



contribute to damage to effector cells, according to the same rate 
as defined by the mathematical model (Table 4). Figure 1(b) shows 
the state chart for the effector cells. In the figure, the cell is either 
alive or dead by age or apoptosis. While the cell is alive, it is also 
able to kill tumour cells and proliferate. In the transition rate 
calculations, the variable TotalTumourCells corresponds to the 
total number of tumour cell agents; and the variable 
TotalEffectorCells is the total number of effector cell agents. In 
the simulation model, apart from the agents, there is also an event 
— namely, treatment - which produces new effector cells with a 
rate defined by the parameter s. 

Experimental design for the simulations. Similarly to the 
experiments from [1], four scenarios are investigated. The 
scenarios have different rates for the death of tumour cells (defined 
by parameter b), effector cells apoptosis (defined by parameter d) 
and different treatments (parameter s). The values for these 
parameters are obtained from [6] (Table 5). In the first three 
scenarios, cancer treatment is considered, while the fourth case 
does not consider any treatment. The simulations for the ABMS 
and the Gillespie algorithm are run five hundred times and the 
mean values are displayed as results. 

Results and discussion. Figure 2 shows the results of our 
experiments. In the first column we display the results from the 
ODE model for guidance. The second column shows the results 
from the Gillespie algorithm and the third column presents the 
ABMS results. Each row of the figure represents a different 
scenario. 

Results for Scenario 1 appear similar for the three approaches, 
although the effector cells curve from the ABMS show more 
variability. To evaluate whether the results are significantly 



different for the two simulation methods, we apply a mixed effect 
model. The nuU hypothesis is that there is no significant difference 
between the methods (and therefore there will be no significant 
fixed effect for the method type). We use a 1 % significance level. 
We are testing the similarity for the population of effector cells. 

The effector cells follow a dynamic similar to 1/x for the time 
between 1 and 100: 



ax(t + b) 



(9) 



We apply a mixed-effect model where the simulation run is 
considered to have a random effect on the parameter a and the 
simulation method has a fixed effect on a and h. The results are 
presented in Table 6. At a 1 % significance level the results of the 
two techniques are significantly different as the p-values for the 
fixed effect of the method on the parameters are less than O.I. 

We believe that the variability observed in the ABMS and 
Gillespie curves, given their stochasticity, also influenced the 
statistical test results. The number of effector cells for all 
simulations foUow a similar pattern, although the similarity 
hypothesis was rejected. This variability of Gillespie and ABMS 
is very evident with regards to the effector cells population as the 
size of the populations involved in the first scenario is relatively 
small, which increases the impacts of stochasticity in the outcomes. 

Results for scenario 2 are shown in the second row of Figure 2. 
The outcomes seem fairly different. By observing the ODE results, 
during about the first ten days, the tumour cells decrease and then 
grow up to a value of about 240 cells, subsequentiy reaching a 
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Figure 3. ABMS state charts for the agents of case 2. 

doi:10.1371/journal.pone.0095150.g003 
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Table 10. Parameter values for case 2. 




Parameter 


Value 


a 


0.18 


b 


0.000000001 


c 


0.05 


aa 1 


g2 


100000 


si 


0 


s2 


0 


mu2 


0.03 


pl 


0.1245 


gi 


20000000 


p2 


5 


g3 


1000 


mu3 


10 


doi:10.1371/journal.pone.0095150.t010 



Steady-State. This initial decrease is also observed in both Gillespie 
and ABMS curves. However, only the Gillespie method shows a 
similar increase in the numbers of tumour cells when compared to 
the ODEs. Similarly to the previous scenario, the Gillespie and 
ABMS simulation curves present an erratic behaviour throughout 
the simulation days. There is, however, an unexpected decay of 
tumour cells over time in the ABMS simulation, which does not 
happen in the Gillespie outcome. We believe the dilference 
observed in the ABMS is due to the the individual characteristics 
of the agents and their growth/death rates attributed to their 
instantiation. While both ODEs and Gillespie are compartmental 
models and therefore they apply the model rates to the cells 
population, ABMS on the contrary, employs these rates in an 
individual basis. As the death rates of the tumour cells agents are 
defined according to the mathematical model, when the tumour 
cell population grows, the newborn tumour cells have higher death 
probabilities, which leads to a considerable number of cells dying 
out. This indicates that the individual behaviour of cells can lead 
to a more chaotic behaviour when compared to the aggregate view 
observed in the compartmental simulation. 

For scenarios 3 and 4, shown in third and fourth rows of 
Figure 2, respectively, the results for the three approaches differ 
completely. The diflFerences are even more evident for the tumour 
cells outcomes. The ODEs results for scenario 3 reveal that 
tumour cells decreased as effector cells increased, following a 
predator-prey trend curve. For the ABMS, however, the number 
of effector cells decreased until a value close to zero was reached, 



while the tumour cells numbers were very different from those in 
the ODEs results. The ODE pattern noticed was possible given its 
continuous character. In the ODE simulation outcome curve for 
the effector cells it is therefore possible to observe, for instance, 
that after sixty days the number of effector cells ranges between 
one and two. These values could not be reflected in the ABMS 
simulation, as it deals with integer values. Similarly, the Gillespie 
approach outcomes did not resemble those from the ODE model. 
There is more variability in the tumour cells curve than in the 
ABMS outcomes, although the number of tumour cells also 
reaches zero after around sixty days. 

In the fourth scenario, although effector cells appear to decay in 
a similar trend for both approaches, the results for tumour cells 
vary largely. In the ODE simulation, the numbers of effector cells 
reached a value close to zero after twenty days and then increased 
to a value smaller than one. For the ABMS simulation, however, 
these cells reached zero and never increased again. For the 
Gillespie model results, a similar pattern as that from the ABMS 
model occurs, although there seems to be less variance in the 
outcome curve. In addition, the mean numbers for tumour cells 
for the Gillespie approach seem smaller that those observed in the 
ABMS. 

Summary. An outcome comparison between an ABMS and 
a Gillespie algorithm model was performed for case study 1. We 
considered an ODE model of tumour cells growth and their 
interactions with general immune effector cells as the baseline for 
results validation. Four scenarios considering small population 
numbers were investigated and results from ABMS and Gillespie 
were different for both populations for all scenarios. Furthermore, 
both approaches differed largely from the original mathematical 
outcomes. These results indicate that, for this case study, the 
stochasticity apphed to the population as a whole when compared 
to that applied to the individual has a higher impact given the 
small population sizes. The result analysis also reveals that 
conceptualizing the stochastic approaches from the mathematical 
equations does not always produce statistically similar outputs. 

Case 2: Interactions between Tumour Cells, Effector Cells 
and Cytokines IL-2 

Case two regards the interactions between tumour cells, effector 
cells and the cytokine IL-2. It extends the previous study as it 
considers IL-2 as molecules mediating the immune response 
towards tumour cells. These molecules interfere in the prolifera- 
tion of effector cells, which occurs proportionally to the number of 
tumour cells in the system. For this case, there are two types of 
treatment, the injection of effector cells or the addition of 
cytokines. 

The mathematical model used in case 2 is obtained from [9]. 
The model's equations described bellow illustrate the non-spatial 



Table 11 


. The results of the mixed model for the 


sequence of local 


maxima in case study 2. 






Technique 


Parameter 


Value 


Std Error 


p-value 


ABMS 


a 


0.122 


0.00073 


0 


Gillespie 


a 


0.131 


0.00102 


0 


ABMS 


b 


442.249 


1.160664 


0 


Gillespie 


b 


432.243 


1.52861 


0 


ABMS 


c 


23149.344 


23.67208 


0 


Gillespie 


c 


22694.685 


31.88635 


0 


doi:10.1371/journal.pone.0095150.t011 
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Table 12. The results of the mixed model for the 


sequence of local minima in 


case study 2. 






Technique 


Parameter 


Value 


Std Error 


p-value 


ABMS 


a 


0.080 


0.00033 


0 


Gillespie 


a 


0.088 


0.00047 


0 


ABMS 


b 


462.004 


1.12224 


0 


Gillespie 


b 


444.888 


1 .46963 


0 


ABIVIS 


c 


17118.133 


25.43450 


0 


Gillespie 


c 


17416.363 


34.44740 


0 


doi:1 0.1 371 /journal.pone.00951 SO.tOl 2 



dynamics between effector cells (E), tumour cells (T) and the 
cytokine IL-2 (//,): 



dE 
dt 



-cT — PL2E 



PiEIl 
gi+h 



(10) 



Equation 10 describes the rate of change for the effector cell 
population E [9]. Effector cells grow based on recruitment {cT) 



and proliferation 



lPlf±L_\ 



The parameter c represents the 



Vi 

antigenicity of the tumour cells (T) [5,9]. ji2 is the death rate of 
the effector cells. p\ and g\ are parameters used to calibrate the 
recruitment of effector cells and s\ is the treatment that wiU boost 
the number of effector cells. 



dT 
~dt 



--a(\-bT)- 



UaET 
' g2 + T 



(11) 



Equation 1 1 describes the changes that occur in the tumour cell 
population T over time. The term a(\—bT) represents the logistic 
growth of T (a and b are parameters that define how the tumour 
cells will grow) and is the number of tumour cells killed by 

effector cells, and g2 are parameters to adjust the model. 



dl, 



PiET 



(12) 



The IL-2 population dynamics is described by Equation 12. 
determines IL-2 production using parameters p2 and ^"3. /.i3 is the 



IL-2 loss. s2 also represents treatment. The treatment is the 
injection of IL-2 in the system. 

Table 7 shows the mathematical model converted into reactions 
for the Gillespie algorithm model. The first column of the table 
displays the original mathematical equation, followed by the 
equivalent reactions and rate laws in the subsequent columns. 

As described in [1], the agents represent the effector cells, 
tumour cells and IL-2. Their behaviours are shown in Table 8. 
The state charts for each agent type are shown in Figme 3. The 
ABMS model rates are the same as those defined in the 
mathematical model and are given in Table 9. In the transition 
rate calculations, the variable TotalTumour corresponds to the 
total number of tumour cell agents, the variable TotalEffector is 
the total number of effector cell agents and TotallLJl is the total 
number of IL-2 agents. In the simulation model, apart from the 
agents, there are also two events: the first event adds effector cell 
agents according to the parameter si and the second one adds IL- 
2 agents according to the parameter s2. 

Experimental design for the simulation. The experiment 
is conducted assuming the same parameters as those of the 
mathematical model (Table 10). For the ABMS and the Gillespie 
algorithm model, the simulation is run five hundred times and the 
average outcome value for these runs is collected. Each run 
simulates a period equivalent to six hundred days, following the 
same time span used for the numerical simulation of the 
mathematical model. 

Results and discussion. The results obtained are shown in 
Figure 4 for tumour cells (left), effector cells (middle) and IL-2 
(right), respectively. As the figure reveals, the results for all 
populations are analogous; the growth and decrease of all 
populations occur at similar times for all approaches. ABMS has 
a little more variability in the results, specially regarding IL-2. We 
believe that for this case, the large population sizes (around lO") 
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Figure 4. Results for the second case study: tumour cells, effector cells and IL-2. 

doi:10.1371/journal.pone.0095150.g004 
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Figure 5. Illustration of the regression models fit for the sequences of the local maxima and local minima for the two different 
simulation techniques. The Gillespie simulations are plotted in purple with the mixed effect models plotted in blue. The ABMS simulation runs are 
plotted in orange with the mixed effect models plotted in red. 
doi:1 0.1 371 /journal.pone.00951 50.g005 



produce a lower variability in the outcomes of the .stochastic 
approaches. For statistical comparison, results were contrasted by 
applying a non-linear mixed effect model, as shown next. 

The local maxima sequence follows a second order polynomial 
fiinction of the form: 



The local minima sequence follows a second order polynomial 
function of the form: 

f(t) = c-a*(t-hf (14) 



f{t) = c + a(t — h)^ (13) For the mixed-effect model, we consider a and h to have fixed 

effects based on the type of simulation (e.g., ABMS or Gillespie 
algorithm) and a and h to have random effects based on the 
individual simulation run. The results of the mixed effect model 



Table 13. Reactions for case 3. 



Name 


Equation 


Rate law (per cell) 


Effector cell recruitment 


-^E: TS 




Effector death 






Effector proliferation 


T + E^E:, T 




Tumour cell growth 


T^lxT 


axT 


Tumour cell death 




K 


Tumour cell death by effector cells 


T + E^E- T 


aaTE 

Si + T 


Tumour growth caused by TGF-/J 


T^2xT;S 


pixSxT 
Si + S 


IL-2 production 


-^I: ETS 


(g4+7-)(l+aS) 


IL-2 decay 






TGV-fi production 




li' + T' 


TGF-^ decay 




ft xS 



doi:l 0.1 371 /journal.pone.00951 SO.tOl 3 
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Table 14. Agents' parameters and behaviours for case 3. 





Agent 


Parameters 


Reactive behaviour 


Proactive behaviour 


Effector Cell 


mill 


Dies 






p\, g\, q\ and ql 




Reproduces 




c 


Is recruited 






cm and ^2 




Kills tumour cells 


Tumour Cell 


a 


Dies 






a 




Proliferates 




aa and g2 


Dies killed by effector cells 






and p2 


Has growth stimulated 






;74 and tetha 




Produces TGF-^ 




c 




Induces effector recruitment 


IL-2 


alpha, p3 and g4 


Is produced 






mill 


Is lost 




TGF-^ 


pA and tetha 


Is produced 






mu3 


Is lost 






p2 and ^3 




Stimulates tumour growth 



doi:10.1371/journal.pone.0095150.t014 



are presented in Tables 1 1 and 12. It can been seen that there is a 
significant difference between the a and b parameter values for the 
two different techniques. We therefore reject the null hypotheses 
and accept that there is a significance difference between the two 
techniques in terms of the sequence of maxima and sequence of 
minima, at a 1 % significance level. Furthermore, the results show 
that the ABMS simulations tend to have larger local maxima and 
smaller local minima, which is clear in Figure 5. 

Summary. Interactions between tumour cells, effector cells 
and the cytokine IL-2 were considered to investigate the potential 
differences and similarities of ABMS and Gillespie algorithm 
outcomes. Statistical comparison between the Gillespie and the 
ABMS results show a significant difference in the outcomes. 
Compared to the original ODE model used as validation, ABMS 
displayed a little more variability in the results, whereas the 
Gillespie algorithm followed mostly the same patterns as those 
produced by the ODEs for all populations in the simulation. As for 
these simulations a bigger number of individuals was required, it 
was also observed that, regarding the use of computational 
resources, ABMS was far more time- and memory-consuming 
than the Gillespie approach. 



Case 3: Interactions between Tumour Cells, Effector Cells, 
IL-2 and TGF-/? 

Case study three comprises interactions between tumour cells 
and immune effector cells, as well as the immvme-stimulatory and 
suppressive cytokines IL-2 and TGF-/? [5]. According to the ODE 
model developed by Arciero et al. in [5], TGF-/? stimulates tumour 
growth and suppresses the immune system by inhibiting the 
activation of effector cells and reducing tumour antigen expres- 
sion. 

The mathematical model is described by the differential 
equations below: 



dE 
~dt 

P\EI 



cT 



Pi 



(15) 



Equation 15 describes the rate of change for the effector cell 
population E. According to [5] , effector cells are assumed to be recruited 
to a tumour site as a direct result of the presence of tumour cells. The 
parameter c in represents the antigenicity of the tumour, 

which measures the ability of the immune system to recognize 
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Figure 6. ABMS state charts for the agents of case 3. 

doi:10.1371/journal.pone.0095150.g006 



PLOS ONE I www.plosone.org 



11 



April 2014 I Volume 9 | Issue 4 | e95150 



Gillespie and Agent-Based for Early-stage Cancer 



Table 15. Transition rates calculations from the mathematical equations for case 3. 



Agent 


Transition 


Mathematical equation 


Transition rate 


Effector Cell 


Reproduce 


p,IE ^ 


p\ X TotallLJ. 








gl + TotallLJ. 



Tumour Cell 



IL-2 
TGF-/ 



doi:10.1371/journal.pone.0095150.t015 



Die 

ProducelL2 
KillTumour 
Reproduce 
Die 

DieKilledByEffector 

ProduceTGF 

EffectorRecruitment 

Loss 
Loss 

stimulates 
TumourGrowth 



p\-- 



c/lS 
ql + Sj 

p3TE 
(g4+T)il+ alphas) 



S2 + T 



T 1- 



^ ] 

1 000000000 y 



r 1-- 



1000000000 



.)) 



..TE 



g2 + T 

p4T^ 
teta^ + 

cT 
1+-/S 

tJ-iS 



plT 



Pl- 



c/1 X TotalTGFBeta\ 
ql+TotalTGFBeta) 



p3. TotalTumour 
{g4+ TotalTumour){\ + alpha.TotalTGF) 

aa X TotalTumour x TotalEffector 
g2 + TotalTumour 

TotalTumour\ 



yTotalTuniour .a yl — 
I TotalTumour. a i 1 — 



1000000000 J 

TotalTumour\ 
1000000000 I 



message from effector 

p4. TumourCells 
teta^ + TumourCells^ 
c 

1 + gamma .TotaltGF 

mu2 

mu3 

pl.TotalTGF 
g3 + TotalTGF 



Table 16. Parameter values for case 3. 



Parameter 



Value 



0.18 
1 



alpha 


0.001 


c 


0.035 


gl 


20000000 


g2 


100000 


g3 


20000000 


g4 


1000 


gamma 


10 


mul 


0.03 


mu2 


10 


mu3 


10 


pi 


0.1245 


P2 


0.27 


P3 


5 


p4 


2.84 


ql 


10 


q2 


0.1121 


theta 


1000000 


k 


10000000000 



tumour cells. The presence of TGF-/J (S) reduces antigen 
expression, thereby limiting the level of recruitment, measured 
by the inhibitory parameter y. The term ji^E represents loss of 
effector cells due to cell death. The proliferation term 

i si+'i ) ip^ ~ 92^) ''^^^^'^'^^ ^^^^ effector cell proliferation depends 
on the presence of the cytokine IL-2 and is decreased when the 
cytokine TGF-j8 is present, pi is the maximum rate of effector cell 
proliferation in the absence of TGF-/J, gi and q2 are half- 
saturation constants, and qi is the maximum rate of anti- 
proliferative effect of TGF-y8. 



dT 



-aTl 1 



OaET pjST 



g2- 



(16) 



Equation 16 describes the dynamics of the tumour cell population. 

The term aT(^l— J) represents a logistic growth dynamics with 

intrinsic growth rate a and carrying capacity K in the absence of 

effector cells and TGF-B. The term ""f^ is the number of tumour 
1^ g2 + T 

cells killed by effector cells. The parameter measures the 

strength of the immune response to tumour cells. The third term 

g^^g accounts for the increased growth of tumour cells in the 

presence of TGF-^. p2 is the maximum rate of increased 
proliferation and is the half-saturation constant, which indicates 
a limited response of tumour cells to this growth-stimulatory 
cytokine [5]. 



doi:l 0.1 371 /journal.pone.00951 50.t01 6 



dl P3ET 
dt^ (g4 + T)i\+aS) 



(17) 
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Figure 7. Results for the third case study: mean values of tumour cells, effector cells and IL-2. 

doi:1 0.1 371 /journal.pone.00951 50.g007 



The kinetics of IL-2 are described in equation 1 7 . The first term 
(g^ + jyi+aS) represents IL-2 production which reaches a maximal 
rate of in the presence of effector ceUs stimulated by their 
interaction with the tumour cells. In the absence of TGF-/J, this is 
a self-limiting process with half-saturation constant g4S [5]. The 
presence of TGF-/J inhibits IL-2 production, where the parameter 
a is a measure of inhibition. Finally, represents the loss of IL-2. 



dS paT^ 



dt + 



(18) 



Equation 18 describes the rate of change of the suppressor 
cytokine, TGF-/J. According to [5], experimental evidence 
suggests that TGF-ji is produced in very small amounts when 
tumours are small enough to receive ample nutrient from the 
surrounding tissue. However, as the tumour population grows 
sufficiently large, tumour cells suffer from a lack of oxygen and 
begin to produce TGF-ji in order to stimulate angiogenesis and to 
evade the immune response once tumour growth resumes. This 

switch in TGF-/J production is modelled by the term p^p^, where 
pn is the maximum rate of TGF-/? production and T is the critical 



tumour cell population in which the switch occurs. The decay rate 
of TGF-/J is represented by the term ^35*. 

Table 13 presents the Gillespie algorithm model used for our 
simulations. The model was obtained by converting the ODEs into 
reaction equations. 

The agents established for the ABMS represent the effector 
cells, tumour cells, IL-2 and TGF-/? populations, as described in 
[1]. The agents' behaviour is defined in Table 14. The state charts 
for each agent type are illustrated in Figure 6. 

The ABMS model rates corresponding to the mathematical 
model are given in Table 15. In the transition rate calculations, the 
variable TotalTumour corresponds to the total number of tumour 
cell agents; the variable TotalEffector is the total number of 
effector cell agents, TotallL^ is the total number of IL-2 agents 
and TotalTGFBeta is the total TGF-/? agents. This model does 
not include events. 

Experimental Design for the Simulation 

The experiment is conducted assuming the same parameters as 
those defined for the mathematical model (Table 16). Similarly to 
the previous case studies, for the ABMS and Gillespie models the 
simulation is run five hundred times and the average outcome 
value for these runs is displayed as result. Each run simulates a 



X 10 



Tumour cells from ABMS 



X 10 
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Figure 8. Results for the third case study: 50 runs for tumour cells. 

doi:1 0.1 371 /journal.pone.00951 50.g008 
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period equivalent to six hundred days, following the time interval 
used for the numerical simulation of the mathematical model. The 
parameters used for the simulations of all approaches are shown in 
Table 16. 

Results and discussion. The mean results of 500 runs for 
the Gillespie algorithm and the ABMS contrasted with the ODE 
model are shown in Figure 7. The left graph in the figure presents 
the outcomes for tumour cells; the graph in the middle shows the 
outputs for effector cells; the graph on the right shows the mean 
IL-2 outcomes (the TGF-/? results have some particularities and 
therefore are discussed next). The figure shows that both Gillespie 
and ABMS do not match properly the original results from the 
mathematical model. Additionally, ABMS is far more dissimilar 
than what was anticipated. In order to understand why the mean 
values were that much difTerent from what was expected, we 
plotted fifty individual runs for each approach, as shown in 



figures 8, 9, 10 and 11. These runs illustrate the variations 
observed in both ABMS (left side of the figures) and Gillespie (right 
side of the figures) approaches, due to its stochastic character. In 
the figures, the ODE model results were also plotted (dashed black 
line) in order to highlight the range of variation produced by the 
stochastic approaches. As it can be observed in the figures, both 
Gillespie and ABMS outcomes produce various slightly distinct 
starting times for the growth of populations. In addition to these 
variations, for a few runs the populations in ABMS decreased to 
zero, as previously reported in [1]. This behaviour was not 
reflected in the Gillespie algorithm results. This indicates that it is 
not always possible to replicate similar results within both 
approaches. 

The use of ABS modelling has therefore led to the discovery of 
additional "rare" patterns, which we would have not been able to 
derive by using analytical methods or the dynamic Monte Carlo 




Days Days 

Figure 10. Results for the third case study: 50 runs for IL2 cells. 

doi:1 0.1 371 /journal.pone.00951 SO.gOl 0 
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Table 17. The results of the mixed model for the 


sequence of local maxima in 


case study 3. 






Technique 


Parameter 


Value 


Std Error 


p-value 


ABMS 


a 


0.1354 


0.001403 


0 


Gillespie 


a 


0.2244 


0.002242 


0 


ABMS 


b 


747.8501 


2.876438 


0 


Gillespie 


b 


595.3515 


3.338526 


0 


doi:l 0.1 371 /journal.pone.00951 SO.tOl 7 



algorithm) and a and b to have random effects based on the 
individual simulation run. The results of the mixed effect model 
are presented in Tables 17 and 18. It can been seen that there is a 
significant difference between the a and b parameter values for the 
two different techniques. We therefore reject the nuU hypotheses 
and accept that there is a significance difference between the two 
techniques in terms of the sequence of maxima and sequence of 
minima, at a 1% significance level. 

Tables 17-18 and Figure 12 show that the sequence of local 
maxima of the ABMS diverge from that of the Gillespie algorithm 
over time. In the ABMS the tumour cells tend to increase to a 
larger count than the Gillespie algorithm simulations causing the 
function of local maxima for the ABMS to be significandy greater 
than the Gillespie. A possible explanation for this, as mentioned 
previously, is the fact that agents have memory and therefore the 
rates (death, proliferation, etc) for a certain cell are determined by 
the cells (and their proportions) present in the system at the 
moment the cell was created. For the Gillespie algorithm, instead, 
the rates are applied globally to the entire population and remain 
constant over the simulation course. Consequently, for Gillespie, 
the individuals do not keep a record of the previous population 
dynamics. This explanation is supported by the observation that 
the function describing the sequence of local minima of the ABMS 
is significandy lower than the Gillespie algorithm over time, as the 
same argument would account for the ABMS simulations reaching 
lower levels. 

Regarding the TGF-/J outcomes, the ODEs results reveal 
numbers smaller than one (Figure 1 1 on the right), which is not 
possible to achieve with the ABMS and the Gillespie algorithm. 
The simulation results regarding these molecules are therefore 
completely different for both stochastic approaches. By observing 
the multiple runs graph of ABMS, however, results indicate that 
the TGF-/? grows at around 100 and 200 days, which resembles 
what occurs in the ODE simulation for the first two peaks of TGF- 
P concentration. This suggests that ABMS, as opposite to 
Gillespie, is capable of capturing some of the behaviours of the 
analytical results even when the outputs are different. We believe 
that these observations need to be further investigated in order to 
determine whether this happens in other case studies. In addition, 
it is necessary to investigate in what circumstances and range of 



Table 18. The results of the mixed model for the sequence of local minima in case study 3. 





Technique 


Parameter 


Value 


Std Error 


p-value 


ABMS 


a 


0.01325 


0.0002059 


0 


Gillespie 


a 


0.01823 


0.0002514 


0 


ABS 


b 


37.33220 


2.1334837 


0 


Gillespie 


b 


5.93249 


2.4013484 


0 
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method, i.e. the Gillespie algorithm. These "extreme cases" found 
by ABMS suggest that there might be circumstances where the 
tumour cells are completely eliminated by the immune system, 
without the need of any cancer therapies. 

We believe that ABMS when compared to Gillespie produces 
extra patterns because of the agents individual behaviour and their 
interactions. While ODEs and the Gillespie algorithm always use 
the same values for the parameters over the entire population 
aggregate, ABMS rates vary with time and number of individuals. 
Each agent is likely to have distinct numbers for their probabilities 
and therefore have its own memory of past events (GUespie, 
however, does not encompass individual memory for its elements). 
The agents individual interactions, which give raise to the overall 
behaviour of the system, are also influenced by the scenario 
determined by the random numbers used. By running the ABMS 
multiple times with different sets of random numbers, the 
outcomes vary according to these sets and the emerging 
interactions of the agents also produce the rare outcome patterns. 

For further statistical comparison of the results that follow the 
same pattern of behaviour for ABMS and Gillespie, a mbced-effect 
model is used. In 236 of the 500 ABMS simulations the tumour 
cell population dies out early in time. The remaining 264 ABMS 
simulations (where the tumour cell population does not die out 
over the [0,600] time period) is compared with the 500 Gillespie 
simulations by a non-linear mixed effect model. 

The local maxima sequence follows a second order polynomial 
function of the form: 

f(t) = Am\\+a{t-bf (19) 

The local minima sequence follows a second order polynomial 
ftmction of the form: 

f{t) = a*{t-bf (20) 

For the mixed-effect model we considered a and b to have fixed 
effects based on the type of simulation (e.g., ABMS or Gillespie 
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Figure 11. Results for the third case study: 50 runs for TGF-/?. The figure on the right is a zoomed version of the figure in the centre. 
doi:10.1371/journal.pone.0095150.g011 



values ABMS is still capable to reflect behaviours of numbers 
smaller than one agent present in the ODE model. 

Summary 

The third case study simulations investigated interactions 
between effector cells, tumour cells and two types of cytokines, 
namely IL-2 and TGF-/?. When compared to the original ODE 
results, both Gillespie algorithm and ABMS produced more 
variability in the outcomes. For each of the five hundred runs, a 
shghdy different start of population growth was observed. In 
addition, ABMS produced extra patterns not observed in the 
original mathematical model and in the Gillespie results. These 
extra patterns have been reported previously in [1] and with the 
present work we wanted to find out whether the Gillespie 
algorithm simulation results would be as informative. This 



indicates that, for this case study, both methods should not be 
employed interchangeably, as some extra possible population 
patterns of behaviour might not be uncovered without ABMS. We 
believe that these emergent examples occur due to the individual 
interactions of the agents and their chaotic character. With these 
results, we answer our third research question that it is not possible 
in this case to obtain extreme patterns using the Gillespie 
algorithm. 

Conclusions 

In this work, we employed three case studies to investigate 
circumstances where we can use ABMS and the Gillespie 
algorithm interchangeably. We aimed at reproducing the 
variability embedded in the ABMS systems to the mathematical 
formulation and verify whether results resemble. Current literature 



o 




0 100 200 300 400 500 

Time 

Figure 12. illustration of the regression models fit for the sequences of the local maxima and local minima for the two different 
simulation techniques. The Gillespie simulations are plotted in purple with the mixed effect models plotted in blue. The ABIVIS simulation runs are 
plotted in orange with the mixed effect models plotted in red. 
doi:10.1371/journal.pone.0095150.g012 
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regarding the comparison of the Gillespie algorithm and ABMS is 
scarce and we wanted therefore to answer the questions: (1) Does 
the Gillespie algorithm produce similar results to ABMS? (2) Can 
these two methods be used interchangeably for our case studies? 
(3) Does the Gillespie algorithm also fmd the extra patterns 
revealed by the ABMS in the third case study? The case studies 
investigated regarded models with different characteristics, such as 
population sizes, modelling effort demanded and model complex- 
ity. 

The first case study involved interactions with general immune 
effector cells and tumour cells. Four different scenarios regarding 
distinct sets of parameters were investigated and in the first three 
scenarios treatment was included. ABMS and Gillespie produced 
different results for all scenarios. It appears that two major 
characteristics of this model influenced the differences obtained: (1) 
The small quantities of individuals considered in the simulations 
(especially regarding the effector population size, which was always 
smaller than ten) that significantly increased the variability of both 
stochastic approaches; and (2) the stochasticity of the Gillespie 
algorithm is applied to the aggregates, while in the ABMS there is 
individual variability. 

Case study 2 referred to the investigation of a scenario 
containing interactions between effector cells, cytokines IL-2 and 
tumour cells. For this case the Gillespie and ABMS approaches 
produced similar outcome curves, which also matched the pattern 
of behaviour of the mathematical model used for validation. As 
populations' sizes had a magnitude of 104 individuals, the erratic 
behaviour of both stochastic approaches was no longer evident in 
the outcomes. However, although results seemed similar, further 
statistical tests reject their similarity hypothesis. It was observed 
that, for case 2 in general, ABMS simulation outcome curves tend 
to have larger local maxima and smaller local minima. 

Case study 3 includes the influence of the cytokine TGF-jS in the 
interactions between effector ceUs, cytokines IL-2 and tumour cells 
from the previous case. The simulation outcomes for the ABMS 
were mostiy following the same pattern as those produced by the 
Gillespie algorithm, although the results were statistically different. 
In addition, Gillespie failed to replicate the alternative outcomes 
found by the ABMS. This indicates that for this case study the 
ABMS results are more informative, as they iUustrate another set 
of possible dynamics to be validated in real-world. Furthermore, 
ABMS was also able to indicate two peaks where TGF-jS 
concentrations have grown, although the corresponding values 
in the mathematical model were smaller than one. 

In response to our research questions, we conclude that 
regarding the interchangeable use of Gillespie and ABMS, 
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